Growth of correlations in gravitational N-body simulations 



Thierry Baertschiger 
Departement de Physique Theorique, Universite de Geneve, 
Quai E. Ansermet 24, CH-1211 Geneve, Switzerland 



o 
o 

(N 
m 



> 
oo 
m 

o ; 
^ . 
o 

6 



Francesco Sylos Labini 
Laboratoire de Physique Theorique, Universite Paris XI, Batiment 211, F-91405 Orsay, Franc^ 

(Dated: February 2, 2008) 

In the gravitational evolution of a cold infinite particle distribution, two-body interactions can 
be predominant at early times: we show that, by treating the simple case of a Poisson particle 
distribution in a static universe as an ensemble of isolated two-body systems, one may capture 
the origin of the first non-linear correlated structures. The developed power-law like behavior 
of the two-point correlation function is then simply related to the functional form of the time 
evolved nearest-neighbor probability distribution, whose time dependence can be computed by using 
Liouville theorem for the gravitational two-body problem. We then show that a similar dynamical 
evolution is also found in a large-scale ordered distribution, which has striking similarities to the 
case of a cosmological CDM simulation which we also consider. 

PACS numbers: 05.25.-a, 45.05.+X, 95.10.Ce, 98.65.-r 



I. INTRODUCTION 

Non-linear gravitational clustering can be studied by 
means of N-body simulations (NBS) which compute nu- 
merically the evolution of a system of particles under the 
action of their mutual gravity. The gravitational many- 
body problem consists in the explanation of the time evo- 
lution of the NBS and in the theoretical understanding 
of the formation of non-linear structures. Up to now, 
two different approaches have been generally studied: on 
the one hand research of approximative solutions of the 
BBGKY hierarchy Q and on the other hand statistical 
thermodynamics mainly developed by Saslaw 2]. 

A main issue in the context of cosmological NBS is 
to relate the formation of non linear structures to the 
specific choice of initial conditions used: this is done in 
order to constraint models with observations of cosmic 
microwave background radiation anisotropics, which are 
related to the initial conditions, and of galaxy structures, 
which give instead the final configuration of strongly clus- 
tered matter. Standard primordial cosmological theoret- 
ical density fields, like the cold dark matter (CDM) case, 
are Gaussian and made of a huge number of microscopic 
mass particles, which are usually treated theoretically as 
a self-gravitating coUisionless fluid H, Q |1| : this means 
that the fluid must be dissipation-less and that two-body 
scattering should be small. The problem then being in 
which limit NBS, based on particle dynamics, are able to 
reproduce the two above conditions. In this context one 
has to consider the issue of the physical role of particle 
fluctuations in the dynamics of NBS as the total energy 
is conserved during time evolution (the only mechanism 



of energy dissipation is related to local gravitational pro- 
cesses). 

In fact, in the discretization of a continuous density 
field one faces two important limitations corresponding 
to the new length scales which are introduced. On the 
one hand a relatively small number of particles are used: 
this introduces a mass scale which is the mass of these 
particles. (In typical cosmological NBS, this mass is of 
the order of a galaxy and hence many orders of mag- 
nitude larger than the microscopic mass of CDM par- 
ticle.) Furthermore, it introduces a new characteristic 
length scale given by the average distance between near- 
est neighbor (NN) particles (A). Clearly the discretiza- 
tion method used should conserve the continuous cor- 
relations, but this is a p roblematic aspects of standard 
methods |E S S S Il0l |- On the other hand one must 
regularize the gravitational force at small scales in order 
to avoid problems related to the divergence of the nu- 
merical integrator and remove coUisional effects due to 
strong scattering between particles. This is usually done 
by using a softening length e in the gravitational poten- 
tial generally defined as 
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This is the second length scale introduced to numerically 
simulate the coUisionless fluid. 

The question which naturally arises is then how to 
choose the two new length scales (A) and e: the flrst 
obvious condition is that they must be both smaller than 
the intrinsic characteristic scales of the continuous fleld 
(as for example smaller than the typical scale correspond- 
ing to the turn-over scale of the CDM power spectrum) . 
Then one has to tune the ratio 77 = (A) /e appropriately 
with respect to the physical problem under study. In 
fact, when rj > 1 one has a larger dynamical range than 
the case ry < 1, but strong scattering between nearby par- 
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tides are not smoothed and hence one is not effectively 
reproducing a dynamics where particles play the role of 
collisionless fluid elements. It is in this sense that one 
talks about the role of discreteness in NBS: that strong 
scattering between nearby particles are produced by the 
discretization and by the choice of 77 > 1 and they should 
be considered artificial and spurious with respect to the 
dynamical evolution of a self-gravitating fluid. This point 
has been considered in different way s and contexts by 
many authors, e.g. [1 S H [H El [H El : they all show 
that discreteness has some influence on the formation of 
the structures. 

For this reason discreteness, which anyway introduces 
large fluctuations in the density fleld up to scales of or- 
der (A), may play an important role in the early stages 
of non-linear structures formation, i.e. when the aver- 
age distance between nearby particles becomes rapidly 
smaller than (A). How discrete effects are then "ex- 
ported" toward large scales, if they are at all, is then 
a deep and difficult problem to be understood. In other 
words the problem is that of understanding whether large 
non-linear structures, which at late times contain many 
particles, are produced solely by the collisionless dynam- 
ics of a fluid and its density fluctuations or whether the 
particle coUisional processes are important also on the 
long-term. For example have argued that discrete- 
ness effects play an important role in the self-similar evo- 
lution of correlated structures, while the effect of NN in- 
teractions has been the subject of a toy model developed 
by 113. 

In [TEIIT^ we have already considered the effects of dis- 
cretization in the dynamics of non linear structure forma- 
tion in several NBS with and without space expansion. 
We have concluded that the fluctuations at the small- 
est scales in these NBS — i.e. those associated with the 
discreteness of the particles — play a central role in the 
dynamics of clustering in the non-linear regime. This 
was based in particular on the fact that the correlations 
appear to be built up from the initial clustering at the 
smallest scales and that the nature of the clustering seems 
to be independent (or at most very weakly dependent) 
on the initial conditions. The theoretical understanding 
of the creation of these correlations should therefore deal 
with the apparently crucial role of the intrinsically highly 
fluctuating initial density field. 

In this paper we put our previous results on a firmer 
physical basis. We study the formation of first structures 
in several NBS. As a reference example we use a cold (zero 
initial velocity) Poisson distribution as initial conditions 
and we consider the case of a non-expanding background, 
i.e. a static universe. In this case, we show that two- 
body interactions are enough to explain the evolution of 
the correlation function at early times, as it has been al- 
ready noticed in E3 • This is done by treating the N-body 
problem as an ensemble of isolated two-body systems. 
Such an approximation is justified, in the Poisson case, 
by the fact that the probability that nearby particles are 
mutually NN is high enough (^^ 0.6) (becoming of order 



one when very close particles are only considered) and by 
the fact that the NN force is the dominating one [l8l| . 
Using Liouville theorem for the gravitational two-body 
problem, we can find the early evolution of the NN prob- 
ability distribution. As this distribution can be linked to 
the conditional density and therefore to the reduced two- 
point correlation function, we also obtain their evolution 
at early times. Comparing with the results from the sim- 
ulations we find an excellent agreement: this shows that 
the first structures observed are a consequence of two- 
body interactions between NNs. After a time of the or- 
der of the typical time scale of two-body interaction, this 
is of course not the case anymore. However we note that 
the functional behavior of the two-point correlation func- 
tion remains unchanged at later times, while the regime 
of strong clustering increases with time. 

We then study in the same perspective three other dif- 
ferent simulations in which the force is not dominated by 
short-scales contributions since the beginning. The link 
between the NN probability distribution is found to be 
an efficient tool to study the nature of the first correla- 
tions developed and the growth of power-law correlations 
when high resolution (77 ^ 1) NBS are considered. 

II. STATISTICAL TOOLS 

A simple tool used to study clustering of a matter 
distribution is the two-point correlation function E3 
(ri(ri)n(r2)) which gives the probability density for find- 
ing one particle around ri and a second one around Y2 
(ri(r) being the microscopic mass density function). In 
the following we will restrict ourselves to distributions 
which have a well-defined average density and are ho- 
mogeneous and isotropic. In that case, the two-point 
correlation function only depends on ri2 = |ri — V2\ and 
the asymptotic average density is positive. This function 
is useful to study both continuous and discrete distribu- 
tions of matter. In the latter case, which is the case of 
interest here, it can be useful to measure averages from a 
point occupied by a particle. For instance, one can define 
the conditional density 

no 

for r > ; this gives the average density at a distance r 
from an occupied point ^ . It is easy to show that one has 
the following relation 

(n(r))p=no[l + e(?-)] forr>0 (3) 

where ^(r) is the non-diagonal part of the reduced two- 
point correlation function jl9j . 



(.)p means that it is a conditional average: the origin is an oc- 
cupied point. 
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In order to study small-scales properties of a discrete 
distribution one may consider the nearest neighbor prob- 
ability distribution uj{r). This gives the probability den- 
sity of the distance from a particle to its NN jlSiJ. Let 
us briefly discuss its relation to the average conditional 
density. By definition, the probability that, given a par- 
ticle, there is another particle in the infinitesimal volume 
element dV^ at distance r is 



pi(r) = {n{r))pdV. 



(4) 



Now we only have to note that the probability uj{r)dr 
for a given particle of having a NN at a distance between 
r and r -f dr is the probability of having no NN in the 
sphere of radius r centered on the particle multiplied by 
the probability of having one particle in the infinitesimal 
spherical shell around this sphere: 

w(r)dr==[l - / uj{s)ds] ■ {n{r))pAT:r^ dr (5) 



where the second part of the right hand side is the prob- 



ability pi (r) with dV = 47rr dr 



where v(r) is given by 



/ ^ 11 3 



(9) 



Averaging on r, we get the probability that two particles 
are mutually NN: 

P3 = / a;(r) exp (— now(r)) dr sa 0.6 . (10) 
Jo 

Hence we have that more than the half of the particles are 
mutually NN. If we restrict ourselves to particles which 
have a NN at a distance / < (A) , this probability becomes 



Pi 



(A) 



uj{r) exp {—nov(r)) dr 



(A) 



(11) 



uj{s) ds 



This result together with the fact that in a Poisson distri- 
bution the force on a particle is mainly due to its NN , 
allows us to consistently treat for an initial short time the 
many-body problem as an ensemble of independent and 
isolated two-body systems. 



III. EVOLUTION OF A POISSON 
DISTRIBUTION 



In the case of a Poisson distribution one simply has 
(n(r))p = no [l9ll. It is then easy to solve Eq. (0) for 
uj{r). One finds^l 



uj{r) = Aimor"^ exp 



(6) 



The average distance between a particle and its NN is 
given by 



(A) 



rLo{r) dr 



ATTUr 



1/3 



(7) 



where F^; is the Euler incomplete gamma function. 

Let us now compute the probability, in a Poisson distri- 
bution, that given a particle and its NN, they are mutu- 
ally NN. Let us suppose that a particle A has the particle 
B as NN at distance r. The probability that A is the NN 
of B is equal to the probability that no other particles are 
in the volume v{r) defined by the portion of the sphere of 
radius r around B which is not contained in the sphere 
of radius r around A. For a Poisson distribution this is 
simply ^ 



P2{r) = exp (-nou(r)) 



(8) 



^ In a Poisson distribution, the probability that there are k parti- 
cles in a volume V is given by (noV)* exp(— noV)/A:!. 



A. Time-scale of NN interaction 

This last result explains what happens if one leaves 
a Poisson distribution without velocity evolving under 
its own gravity: most of particles will fall on their NN. 
Let us determine the time-scale of this phenomenon. To 
this aim, one can use conservation of energy in a pair of 
particles of mass m: 



E = -- 



Gr 
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where we have used the Newtonian potential. As we will 
see in more detail in the next subsection, the problem can 
be reduced to a single dimension and choosing center of 
mass coordinates, we get xi{t) = —X2{t). After some 
algebraic manipulations Eq. H12() becomes 



(13) 




assuming that 2:1(0) > 0. The time of fall is 

-1/2 



tfi,ii{ro) = - 
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(14) 



Taking for rg the mean distance between NNs, (A) given 
by Eq. Q, we get 



^=T\/3r|(4/3). 



1 



1.148 



where po = TnriQ is the mass density. 



(15) 
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B. Approximate evolution of the conditional 
average density 

As already mentioned, the force on a particle in a Pois- 
son distribution is almost only due to its NN. As in our 
simulations the particles have no initial velocity, they will 
start to fall in direction of their NN and we will see that 
this is what explain the early evolution of {n{r))p for a 
time t ^ T. 

Let us consider that the interaction potential is U (r) = 
U{r) As said before, we make the assumption that the 
force on a particle is only due to its NN and that the 
Poisson distribution can be approximated by an ensemble 
of particle pairs evolving independently. The evolution 
of one of these pairs is given by the following equations: 



mxi = -VxiJ/(ri2) 



dU 




Xl - X2 


dr 






dU 




Xl - X2 


dr 


ri2 


ri2 



(16a) 
(16b) 



with ri2 = |xi — X2I. Adding these two equations, one 
gets Xl = — X2 (conservation of total momentum). As 
the particles are supposed to be at rest at i = 0, one has 
Xl = — X2 with a proper choice of the origin (center of 
mass coordinates). With this relation, xi — X2 = 2xi = 
—2x2 and one has to solve only one equation of motion, 
for particle 1 for instance: 



mxi 



dU 
dr 



2|xi| 



Xl 

|xi| 



(17) 



Using again the fact that the initial velocity is null, one 
can reduce the number of dimensions to one: 



dU 
dr 



2\x\ 



dV 

sign(x)=-— (18) 
dx 



with V{x) = U{\2x\)/2 which is the equation for the 
evolution of a single particle in the potential V. One 
can now use Liouville theorem in order to study the 
evolution of a phase space density function of systems 
evolving according to this equation and choosing an ap- 
propriate density function, one can obtain the evolution 
of uj{r). 

If /(x, ti, t) is a phase space density function, Liouville 
theorem states that 



[dt + vd, + id,)f = 



(19) 



where, in our case, mv — —dV/dx. The appropriate 
initial condition is 



(20) 



^ We do not restrict ourselves to a precise potential as it can vary 
in different NBS. 



with ijj{r) given by ©. We divide it by 2 in order to 
have half of the particles with a; > and half with x < 
0. Knowing f{x,v,t), one can obtain Lu{r,t), the time 
evolved NN probability distribution, with 

/OC POO 
f{-r/2,v,t)dv+ / f{r/2,v,t)dv. 
-00 'J —00 

(21) 

In order to solve the Liouville equation, let us denote 
by (t)t{xo,vo) = {Xt{xo,vo),Vt(t,xo,vo)) the solution of 
the equation of motion with initial condition xq, voatt — 
0. The Liouville equation implies that f{x, v, t) remains 
constant along a phase space trajectory: 



} {Xt{xQ,vo),Vt{xQ,va),t) = /(xcWcO). 



(22) 



With our initial conditions, the solution of this equation 
is therefore 

f{x,v,t) = /((/)_f(a;,u),0) 

dxidwi [/(0_f(a;i,i;i),O) 

X S{x — Xl) S{v — vi)] 

dxoduo [f{xo,VQ,0) 

X S{x - Xt{xo,vo)) S{v - Vt{xo,vo))] 

(23) 

with f{Xt{x,v),Vt{x,v),t) = f{(t)t{x,v),t). We have 
used the fact that the determinant of the Jacobian ma- 
trix in the change of variables from (xi^vi) to (xo,wo) 
is det{d(l)t/d{x,v)) = 1. This is actually related to the 
Liouville theorem |2ll |. 

With this solution, we can get the evolution of u!{r,t). 
First let us remark that f{x,v,0) = f{—x,—v,0) and as 
the force in (|18|l is odd, if x{t) is a solution, —x{t) is 
also a solution. This permits to show that f{x,v,t) = 
f{—x,—v,t). It is then easy to see that (PT|) can be 
rewritten as 



io{r,t)^2 f{r/2,v,t)dv. (24) 
>.' — 00 

Using this last equation and Eq. H23|l . one has: 

w(r,t) = 2 



dxo / duo /(a;o,«o,0) 

x5(r/2-Xt(xo,«o))- (25) 
As f{x,v,0) — S{v)uj{2\x\)/2, this becomes 

POO 

Lo{r,t)= I dxouo{2\xo\)5{r/2-Xt{xo,Q)). (26) 

one has 

(27) 



Using the fact that for a function / : R — 

5{x - y) 



y&z{f} 



\f'iy)\ 
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with Z{f) = {y£ 
Lo{r,t) = 



E 



f{y) = 0}, we get 



u;{2\xo\) (28) 



with — {xq £ R I Xt{xo,0) — r/2}. Of course there are 
points xo in SI such that dXt {xq , 0) /dxo = and there- 
fore w(r, t) is not weU defined at some isolated points. 

We may solve numerically Eq. (|18() to find a solution 
for Xt{xQ, 0). The steps to get uj{r, t) for a given t are the 
following. We start with a set {xo,i — a^o.min -\- i5 \ 5 > 
0, < z < rt} where Xmim ^ and n are chosen so that the 
region covered gives non-negligible values for lo{2x) and 
that this region is sufficiently sampled. For each i, one 
calculates numerically Xi = Xt{xQ,i^ 0). By doing a linear 
interpolation with these values, we have an estimate of 
Xt{xo, 0) for all Xq in the region covered by the xq^. The 
last step is to find the set of x which solve Xt{x,0) = r/2. 
Once we have w(r, t), we find the conditional density by 
using Eq. ©. 
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FIG. 1: The normalized conditional density in a Poisson dis- 
tribution at time r, 2r, 3r, 4r. Note that once correlations 
are developed, the subsequent evolution increases the range of 
scales where non-linear ((n(r))p 3> no) clustering is formed, 
while the function behavior of two-points remain unchanged. 



C. Comparison with a simulation 

In order to test the simple argument presented in 
the last subsection, we did a N-body simulation. We 
have used the code Gadget based on a tree algo- 
rithm. The infinite universe is simulated by using pe- 
riodic boundary conditions and the usual Ewald sum- 
mation technique. The force between two particles is 
not exactly Newtonian but a softened one is used [2^ . 
Note that the potential used is not the standard Plummer 
one but a similar one which has the advantage of being 
perfectly Newtonian at a scale larger than the softening 
length. 

We have generated a Poisson distribution with N ~ 
2)2^ particles in a box of nominal size L. The mass of 
the particles is such that the mass density is one. The 
softening length is e = 0.00175L: by using Eq. Q we find 
(A) « 0.017L and hence rj « 10. The initial velocities are 
set to zero, and the simulation is run up to 4 r. 

The time evolution of the conditional density is shown 
in Fig. ^ (here and in what follows we normalize the con- 
ditional average density to the asymptotic density, i.e. 
we consider {n{r))p/no). It is worth noticing that once 
the power-law correlations are developed, the subsequent 
evolution increases the range of scales where non-linear 
clustering is formed, i.e. where {n{r))p ^ hq, by approx- 
imatively a simple rescaling: denoting by {n{r, t))p the 
conditional density at time t, one has 

{n{r,t + 5))pK {n{a-r,t))p (29) 

where a > is a constant which depends on the time |l6j | . 

In Fig. 121 we show the initial NN density distribution 
obtained from the Poisson distribution used in the simu- 
lation and the one from Eq. © . The conditional density 
of the initial configuration together with the one obtained 
by using Eq. (|SJ| are shown in Fig. |31 
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FIG. 2: Initial NN probability distribution for the Poisson 
case. The solid line is the exact solution given by Eq. © 
while the dashed line is the measured one in the simulation. 

The evolution of the NN probability distribution in 
the simulation together with the one obtained from H28|l . 
at times 0.5 r, r and 1.5 r, is shown in Figs 14161 We 
may notice that the agreement is quite good and even 
excellent at t = 0.5 r. The differences which appear at 
t = T and t — 1.5 r seem to be explained by the following 
arguments. 

First of all we remind that in a Poisson distribution 
the force acting on a particle can be decomposed in two 
terms: the one given by the NN particle and the one 
due to all the other particles. While the first represents 
a large contribution, the second rapidly goes to zero for 
symmetry reasons [l3 |. However, for particles which have 
a NN further than the average (A), the situation is dif- 
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FIG. 3: Normalized two-point conditional density at the ini- 
tial time: the solid line (Theory) is the theoretical ensemble 
average behavior while the dashed line is the measured one in 
the simulation. 



ferent. Let us denote by A such a particle and its NN 
by B. The force Fba from the latter on A being weaker 
than the average force on a particle from its NN, the 
force contribution of other particles nearby becomes also 
important on the total force on A. This total force is 
then not necessarily in the direction of B and the parti- 
cle A will not "fall" on it. Furthermore the particle B 
has a high probability of having a NN different from A] it 
should therefore not go towards A. The distance between 
A and its NN B should then grow. This is actually what 
we observe if we compare carefully Figs El and El looking 
the value of r/L at large scales at which the NN prob- 
ability distribution reach a value of 10~^, we see that it 
is 3 • 10-2 at i = and 3.5 • lO'^ at t = 1.5 r, i.e. the 
particles whose NN is at a distance 3 • lO^^L initially are 
at a larger distance (3.5 • 1Q~^L) at i = 1.5 r. 

Secondly, concerning particles which have their NN at 
a distance closer than the average (A) we observe that at 
scales between \Q~^L and 5 • \Q~^L a bump is created: 
our simple model predicts less particles than observed in 
the simulation. This seems to be a sign of the creation 
of larger structures. If two particles are isolated, they 
will move in a regular oscillating motion. This is what 
the model predicts. In the simulation these two particles, 
i.e. a particle and its NN, will move together for a while 
as in the model but in the same time be attracted toward 
another pair or group of particles, which is not described 
by the model. This could have the effect of bringing 
the two particles closer together and even give rise to 
an exchange of NN with the other group of particles, 
making the evolution of the NN probability distribution 
evolving differently from the model. The bump reflects 
therefore this step of the clustering which tends to bring 
pairs together. 

In Figs I7I9I we compare the predictions of the condi- 
tional density given by our model, with the measured 



FIG. 4: NN probability distribution aXt = 0.5 t in the Poisson 
simulation. 
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FIG. 5: NN probability distribution at f = r in the Poisson 
simulation. 
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FIG. 6: NN probability distribution at f = 1.5 r in the Poisson 
simulation. 
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FIG. 7: Behavior of the average conditional density at t = 
0.5 r in the Poisson simulation. 

ones from the simulations at times 0.5 r, t and 1.5 r. 
One sees that our approximation works again really well 
as it succeeds in reproducing the development of the cor- 
relations. This means that these correlations are there- 
fore only a consequence of the interaction of NNs. We 
may also notice an interesting thing: aX t = 1.5 r, even 
if the agreement is marginally good at scales larger than 
10~^L, it is still correct at smaller scales. An explanation 
is that these scales correspond to pairs whose particles 
were very close (i.e. < (A)) initially and therefore well 
bounded. When they start to feel the effect of particles 
around, their relative motion is not affected and is still 
described by a two-body interaction. 

At larger scales, where there are no correlations, our 
approximation fails to reproduce the correct behavior at 
all times. For a certain r the conditional density goes 
rapidly to 0. This is due to the fact that at these scales, 
the NN probability distribution is really small and the 
Eq. ISJ is not valid anymore: this equation implies that 
the density around a particle is only due to its NN and 
that there are no particles further than the NN. There- 
fore at distances larger than the average distance between 
NNs, the density has to go to as there are no other par- 
ticles to maintain a non-zero density. 

We finally remark that, as noticed in to study the 
role of these NN interactions in the evolution of cluster- 
ing, one may modify the force integrator of the numerical 
code to include only the NN contribution to the gravita- 
tional force. Of course, the result agrees perfectly with 
the study presented here. 



D. Poisson with large softening 

In the Poisson simulation, we have observed that the 
first structures created are pairs of particles. Now we 
present another simulation in which this is not the case. 
It is simply a Poisson simulation with a large softening. 
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FIG. 8: Behavior of the average conditional density aX t = t 
in the Poisson simulation. 
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FIG. 9: Behavior of the average conditional density at t = 
1.5 r in the Poisson simulation. 



one hundred times larger than in the previous case: e = 
0.175L and hence 77 ~ 0.1. 

Figure [TUl shows the evolution of the conditional den- 
sity in this simulation. The time is still in unit of r but 
only for comparison with the first Poisson simulation be- 
cause this is not anymore a microscopic characteristic 
time. One can see that the correlations do not develop 
at the smallest scales of the system ((A) = 0.017i) but 
are directly found up to 10~^L which is of the order of e. 

Looking now at Figs II 11141 where we compare the 
conditional densities obtained from the simulation and 
the ones reconstructed from the NN probability distribu- 
tions, we see that as soon as correlations develop they 
are already made of many particles as the approxima- 
tions of the conditional density by the NN probability 
distribution fails. 

In the beginning of this simulation, the NN contribu- 
tion to the total force acting on a particle is clearly not 
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FIG. 10: Evolution of the normalized conditional density in 
the Poisson with large softening simulation. The times are 
0, 1, 2, 3, 4 in units of t. 
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FIG. 11: Reconstruction of the conditional density from w(r) 
at t = in the Poisson with large softening simulation. 



important. The dominant contribution is actually the 

force due to many particles at some larger scales. This 
means that two nearby particles do not fall on each other 
as in the previous case but feel approximatively the same 
force and therefore go in the same direction once the sim- 
ulation starts. This direction should be the one of the 
nearest mass over-density. Some other particles will also 
be attracted in this direction. The effect of this motion 
is the formation of the first structures, directly made of 
more than two particles. 

As a final remark, it is interesting to note that when 
power-law correlation are formed at t « 4t the exponent 
and the amplitude of the conditional density agrees very 
well with the simulation with small softening previously 
discussed. 



FIG. 12: Reconstruction of the conditional density from a;(r) 
at i = r in the Poisson with large softening simulation. 
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FIG. 13: Reconstruction of the conditional density from u){r) 
at i = 2 r in the Poisson with large softening simulation. 
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FIG. 14: Reconstruction of the conditional density from a;(r) 
at t = 3 r in the Poisson with large softening simulation. 
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IV. THE SHUFFLED LATTICE AND THE CDM 
CASE 

We study now two different cases where the average 
force on a particle in the initial distribution is different 
from the Poisson case, i.e. it is not dominated by the 
NN one. The main point here is to use the relation © 
to study the creation of the first structures: by obtain- 
ing the NN probability distribution in a simulation, we 
reconstruct the conditional density and compare it with 
the one measured directly in the simulation. 



A. The shuffled lattice 

A shuffled lattice is a simple ordered distribution 
which is obtained by adding a random small perturba- 
tion to a perfect lattice of particles: each particle of this 
lattice is moved randomly in a cubic box centered on the 
unperturbed position of the particle. The only parameter 
is then the ratio 



(30) 



between the size of the cubic box 26 and the lattice spac- 
ing I. When as = 0, it is a perfect lattice while as Os ^ cxd 
it becomes a Poisson distribution 19]. For the simula- 
tion presented here, we have used a shuffled lattice with 
32'^ particles, and shuffling parameter — 0.25. The 
mass of the particles, the number density and the soft- 
ening length of the force are the same as for the Poisson 
simulations previously discussed: this gives 77 « 14. 

In Fig. 1151 the evolution of the conditional density is 
shown. The time goes from to 4 t with r given by 
Eq. (|15|l . One may note that once correlations are devel- 
oped, the evolution proceeds in a very similar way to the 
Poisson case (see Fig. ^| : it is the same kind of rescaling 
as given by Eq. H29|) . the only difference being the speed 
at which this happens 16] . 

In Fig. ^1 it is shown the NN probability distribution 
measured in the simulation at the corresponding times. 
It is important to note that at t = there is an anti- 
correlation at small scales: the normalized conditional 
density is smaller than 1. This is due to the fact that 
two particles cannot be closer than a minimal distance 
which depends on as', the excluded volume feature be- 
ing a typical property of super-homogeneous distributions 
[l9l ] . This can be seen with the NN probability distribu- 
tion which is very peaked around the mean inter-particle 
separation. 

The Figs 1171 to 1211 show the reconstructed conditional 
density (by using the NN probability distribution) and 
the one measured directly in the simulation. As for the 
Poisson case, one sees that the first structures observed 
via the conditional density are only due to a change of the 
NN probability distribution. Of course the dynamics of 
a particle with its NN are not described in the same way 
as for the Poisson case. The force on a particle cannot 
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FIG. 15: Evolution of the normalized conditional density in 
the shuffled lattice distribution. The times are 0, 1, 2, 3, 4 in 
units of r. 
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FIG. 16: Evolution of the NN probability distribution in the 
shuffled lattice for the same times than in Fig. 1151 



be approximated by the one from its NN but the latter 
seems to be sufficiently important to give the direction 
of the particle displacement. Another interesting point is 
the fact that there are two phases in the clustering. This 
can be seen in Fig. E| between tQ — and ii = t almost 
nothing happens while between ti — t and t2 = 2t the 
correlations are quickly developed. As t2 — ti — t is the 
typical time-scale for two isolated particles, separated by 
a distance of order I, to fall on each other, this seems to 
show that this brutal change is a sign of such a behavior. 

In order to verify this argument we have done a sim- 
ple test: we have run the simulation again but with a 
modified integrator which, for a given particle, calculates 
the force acting on it only from its n NNs, n being an 
integer identical for all the particles, which can be chosen 
arbitrarily and changed during the simulation. For our 
study what we have done exactly is the following: 
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FIG. 17: Reconstruction of the conditional density from u>{r) FIG. 20: Reconstruction of the conditional density from uj{r) 
at t = in the shufHed lattice. at t = 3 r in the shuffled lattice. 
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FIG. 18: Reconstruction of the conditional density from uj{r) 
at t = T in the shufHed lattice. 



FIG. 21: Reconstruction of the conditional density from uj{r) 
at t = 4r in the shuffled lattice. 
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FIG. 19: Reconstruction of the conditional density from Lu{r) 
at t = 2 r in the shuffled lattice. 



1. at t = 0, the integrator finds the 6 NNs of each 
particle; 

2. it starts to evolve the system up to t = r but at 
each time step, the force on a particle is due only 
to the 6 particles, which were its 6 NNs a.t t — 0; 

3. at t = T, the integrator finds the NN neighbor of 
each particle; 

4. it continues the evolution up to i = 2 t, the force on 
a particle being now only the one from the particle 
which was its NN a.t t — t. 

In Fig. 1221 we show the result which confirms our as- 
sumption: between r and 2 r, the dynamics is driven by 
NN interaction. Furthermore one can see that between 
and T, what matters for a particle is the force from its 
6 NNs chosen for the reason that in a perfect lattice, for 
a given particle, there is not a single NN but there are 
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FIG. 22: Conditional density in the simulation (S) and in the 
modified simulation (MS) for the shufiled lattice. 



6 NNs, all at the same distance (A). In the case of the 
lattice, these 6 particles are in a perfect symmetric con- 
figuration around the center particle (this is the case for 
all the particles when considered as centers) . This implies 
that the force resulting from these particles cancels. In 
a shuffled lattice, as long as the parameter as is smaller 
than one, this remains the case: even if there is a single 
NN for each particle, there are always 5 others particles 
which are almost at the same distance as the NN. The 
simple test presented shows actually that the force from 
these 6 particles is what matters for the evolution of the 
correlations in the system between t = and r and the 
force from more distant particles is negligible. 

Some simple calculations show actually that the force 
on a particle in a shuffled lattice is approximatively given 

by 



F3 = 2^/3^ 



(31) 



assuming Gm^ — 1. Looking at Fig. 1231 one has for in- 
stance that the squared distance vqi between the central 
particle and particle 1 is given by 



'01 



E 

k=y,z 



(32) 

where Si = {si^x-, £i,y, £i.z) is the displacement of the zth 
particle with respect to its lattice point. Supposing that 
these displacements are small compared to I, one finds 
that the force on the central particle from particle 1 is 



Fi.k 



(33a) 



+ 0{e^) iovk^y,z, (33b) 



with Ei^k = £i,kll- Making now the sum over the 6 
particles around and averaging on all the e^.^ which 



are random variables going from —5 to 5, one obtains 
{Fx) = {Y![Fi,x) = and a variance {F"^) = 46'^ /l^. 
This gives for the total squared force 



{F') = {F-^ + F^ 



Ft 



12^ 



(34) 



whose square root is given by Eq. (|31|l . The force from 
the NN is given approximatively by -Fnn ~ which 
shows that the real force is roughly 



w 273 asFi 



(35) 



One can estimate a time scale ts defined by the relation 
1/2 — GmFst1/2 which is an approximative upper bound 
for the time scale needed by two NN particles to fall on 
each other: 



1 



2ti 



V3 flg ^/i7^Gpo 



1.7- 



(36) 



with T given by Eq. IjlSfl . Some simple numerical tests 
performed by varying show that the real time scale is 
closer to 



(37) 



In our case, with Og = 0.25, this gives Ts « 2r which is 
indeed in good agreement with the simulation. 
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FIG. 23: Force in a shufiled lattice. The small circles (o) 
represent the lattice points while the black dots (•) represent 
the particles. 

In summary, while in a Poisson simulation, the corre- 
lations are made directly from the interaction between 
NN, in a shuffled lattice, there is a first phase in which a 
particle interacts mainly with his 6 NNs. This phase is 
characterized by strong anti-correlations which are slowly 
destroyed. This is then followed by a second phase in 
which some positive correlations are rapidly developed 
under some dynamics driven by NN interactions. 



B. CDM simulation 

Finally, we study a CDM cosmological simulation 
which has been done by the Virgo Consortium 0| . This 
simulation is representative of many other cosmological 
simulations as their parameters, their initial particle con- 
figurations and their small scales properties are more or 
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less always the same. The following discussion should 
therefore apply to other cosmological simulations of CDM 
type. Compared to the simulations we have done, this 
cosmological simulation is different on two points. Firstly 
there is space expansion. Secondly the initial conditions 
(IC) are very elaborated. This last point needs some ex- 
planations. 

The goal of this simulation is to study the evolution of 
a gravitating fluid made of CDM particles with particular 
initial correlations. As already mentioned, the particles 
in the simulation do not represent CDM particles but are 
kind of clouds of CDM whose mass are of the order of a 
galaxy. This discretization of the fluid introduces some 
effects which are reduced by putting initially the particles 
in a particular way. The trick is to create first a distribu- 
tion where the force on a particle is almost zero. In the 
Virgo case, this is done by running the integrator used 
for the simulation on a Poisson distribution with a neg- 
ative gravity constant during a while. The distribution 
obtained is characterized by the fact that the main part 
of the force on a particle comes from large scales mass 
fluctuations. The contribution from nearby particles is 
negligible. Note that the use of a repulsive gravity gives 
a behavior similar to a one component plasma Q . 

On this new distribution, it is necessary to apply a 
correlated displacement which would transform a contin- 
uous and perfectly uniform distribution into the expected 
CDM fluid, i.e. with a power spectrum on relatively large 
scales ^ behaving as P{k) ^ /c" with —3 < n < —1. This 
displacement field is applied by using the Zeldovich ap- 
proximation which also fixes the initial velocity of each 
particle as a function of its displacement. The distribu- 
tion obtained is therefore correlated at all scales and has 
some small initial velocity ^ . 

Note that as the pre-initial distribution has super- 
homogeneous properties (as a lattice or a one-component 
plasma 0, there are two main points to be con- 
sidered: (i) on small scales the distribution continues 
to have the excluded volume feature typical of super- 
homogeneous systems, {ii) on large scales the correlations 
properties are given by a complex combination of the pre- 
initial correlations (which are long-ranged) and by the 
correlations imposed by the displacement field. Whether 
the resulting fluctuations field has the same small-scales 
properties of the CDM continuous distribution is ques- 
tionable 0, 0, H, E|. However here we are interested 
only on the small-scales features which have the clear 
imprint of the pre-initial super-homogeneous distribution 
very similar, as we discuss below, to the shuffled lattice. 

This simulation is made with N = 256^ particles in a 
box of size L = 239.5 Mpc//i (where 0.5 < ft, < 1 is the 



* It is not the aim of this simulation to consider the small k region 
where P{k) ~ k. 

^ Note that the velocities are small I3| . This is why afterwards we 
dare to compare this simulation with our initially static simula- 
tions. 
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FIG. 24: Evolution of the NN probability distribution in the 
CDM simulation. The times are given by the redshift z. 
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FIG. 25: Evolution of the conditional density in the CDM 
simulation. The times are given by the redshift z. 



dimensionless Hubble constant). The masses are such 
that = 1 and it should represent a standard CDM 
model. The softening is e = 0.036 Mpc//i which gives 
•q w 25. This simulation goes from a redshift z = 50 to 
z = 0. 

We have measured the conditional density and the 
NN probability distribution. With the latter we have 
used the approximation based on NN probability distri- 
bution to compute the conditional density. The results 
are shown in Figs 1261281 while the evolution of the NN 
probability distribution and the conditional density are 
shown in Figs |^ and 

The first striking feature that we note is that the evolu- 
tion is very similar to the shuffled lattice case. The con- 
ditional density from being anti-correlated distribution 
develops positive power-law like correlations at scales 
smaller than (A). 

This evolution of the correlations is well described by 
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FIG. 26: Reconstruction of the conditional density from u>{r) 
in the CDM simulation ai z — 10. 




r/L 

FIG. 27: Reconstruction of the conditional density from Lu{r) 
in the CDM simulation at z = 5. Note that the discrepancy 
at scales below e comes from a too small statistics on the 
measured conditional density. 



using the NN probability distribution, which means that 
these correlations are simply due to correlations between 
NNs. In we have already analyzed this simulation. 
We had observed that correlations started at the smallest 
scales of the system, i.e. e < r < (A). Now with the 
relation between the NN probability distribution and the 
conditional density, we can make this observation more 
accurate: the "correlations at the smallest scales" are 
actually correlations between NNs. As in we can 
again raise the question of whether these correlations are 
due to some interactions between NNs or are a "large 
scale" effect, i.e. a consequence of the initial velocity field 
and the acceleration of the particles under the gravity 
of large scales mass fluctuations. This large scale effect 
would be what we expect from fluid dynamics. 

The main point in Tsj and (l^ was the kind of uni- 
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FIG. 28: Reconstruction of the conditional density from uj{r) 
in the CDM simulation at z = 3. Note that the discrepancy 
at scales below e comes from a too small statistics on the 
measured conditional density. 

versality of the correlations developed in different gravi- 
tating systems of particles, among them this CDM simu- 
lation, a Poisson and a shuffled lattice. Now we can add 
that the first correlations are exactly of the same kind in 
all these simulations, namely NN correlations. As in a 
Poisson and a shuffled lattice, the discretization plays an 
important role in the creation of these correlations, this 
would suggest that it is the case for the CDM simulation. 



V. CONCLUSIONS 

The fundamental relation used in this paper is Eq. |(SJ). 
It relates the NN probability distribution aj(r) to the con- 
ditional density (n{r))p at scales of the order of the aver- 
age distance between NNs as long as most of the particles 
have a clear NN. By checking if this relation holds in a 
simulation, we get an interesting information on the na- 
ture of the correlations: are they only due to NN corre- 
lations or do they show the existence of structures made 
of many particles. 

In three simulations that we have considered, Poisson, 
shuffled lattice and CDM, which are high-resolution ones 
(77 S> 1), we have seen that this relation holds at early 
times showing that the correlations grow by being ini- 
tially only NN correlations. In another simulation, Pois- 
son with large softening such that ^ 1, we have seen 
that this is not the case anymore. In this simulation, 
the first correlations are due directly to the formations 
of large structures — i.e. larger than the typical distance 
between NNs — containing more than two particles. 

The results for the high-resolution Poisson simulation 
and for the shuffled lattice has encouraged us to push 
the analysis a bit further. In the Poisson case, using the 
relation ||SJ) and considering the following facts: (i) the 
force on a particle is mainly due to its NN, (ii) for more 
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than half of the particles, two particles are mutually NN; 
we could treat the system as a set of isolated two-body 
systems. Knowing the initial NN probability distribution 
and using the Liouville theorem, it has been then possible 
to find the early evolution of the correlations with quite 
a lot of precision. 

In the shuffled lattice, we have observed that at the 
beginning the situation was more complex than in the 
Poisson case. Due to the approximative symmetries, the 
early evolution involves interactions between more than 
two particles. Actually, in this case, instead of having a 
single NN for each particle, there are six particles which 
lie at almost the same distance: such a situation changes 
the small scale behavior of the force on an average par- 
ticle by introducing a compensation, which is exact and 
gives a null force only when the lattice case is consid- 
ered. However, the result is similar to the Poisson case: 
the formation of correlations between NNs. At a later 
time, the situation becomes exactly identical to the Pois- 
son case as the system behaves like a set of isolated two- 
body systems. The consequence is a rapid growth of the 
correlations between NNs. 

For the CDM simulation, we have not tested whether 
the evolution could be explained at a certain time by NN 
interactions. This is a really important question because 
this simulation is supposed to describe a fluid. The parti- 
cles arc not meant to describe particles but mass tracers: 
they should follow the flow due the gravity from the large 
scales mass fluctuations. If this simulation would have 
the same dynamics as the Poisson and shuffled lattice, 
that is that it could be explained by NN interactions dur- 
ing a small amount of time, this would clearly show that 
the fluid is not well simulated as the evolution would be 
influenced by the discrete nature of these particles result- 
ing from the discretization of the fluid and which would 
therefore not exist in a real fluid. This would then re- 
quires some careful studies in order to understand how 
these effects influence the later evolution. 

In some previous papers 0,^^, we had already raised 
these questions, after having observed the kind of uni- 
versal correlations developed in different simulations all 
characterized by their particle based dynamics. In some 
recent papers 0, j some others authors have also an- 
alyzed the influence that these particles could have on 
the evolution but on the consequences of close encoun- 
ters between these particles. Their conclusions were that 
it has an influence on the density profile of the clusters. 

With this paper we have tried to bring a new element 
to the understanding of what happen in such several high 
resolution simulations, including the cosmological CDM 
one, by showing the nature of the first correlations devel- 
oped but we also raise some new questions which should 
clearly deserve further studies. From our results we now 
argue for three conclusions about the nature of cluster- 
ing in the non- linear regime observed in these NBS. With 
respect to cosmological NBS, we conclude that the ex- 



ponent characterizing the non-linear clustering observed 
has essentially nothing to do with (z) the expansion of the 
Universe, or (ii) the nature of the small initial fluctua- 
tions imposed in the IC. We further present evidence for 
the qualitative description of the dynamics driving this 
clustering given in [TJ based on the Poisson case, and in 
101 based on a similar analysis of the CDM simulation: 
(in) the non-linear clustering develops from the large 
fluctuations intrinsic to the particle distribution at small 
scales (specifically around the smallest resolved scale e). 
In particular we show here that the exponent charac- 
terizing it can be seen to emerge at early times in the 
simulations when the evolution is well approximated as 
being due only to the interactions between NN particles. 

A more quantitative description of this dynamics is ev- 
idently needed, with the principal goal of understanding 
the specific value observed of the exponent. In the cosmo- 
logical literature (see e.g. 0) the idea is widely dispersed 
that the exponents in non-linear clustering are related to 
that of the initial power-spectrum of the small fluctua- 
tions in the CDM fluid, and even that the non-linear two- 
point correlation can be considered an analytic function 
of the initial two-point correlations HE Ell (although, see 
where more emphasis is put on the tendency for IC 
to be washed out in the non-linear regime). The mod- 
els used to explain the behavior in the non-linear regime 
usually involve both the expansion of the Universe, and a 
description of the clustering in terms of the evolution of 
a continuous fluid. We have argued that the exponent is 
universal in a very wide sense, being common to the non- 
linear clustering observed in the non-expanding case. It 
would appear that the framework for understanding the 
non-linear clustering must be one in which discreteness 
(and hence intrinsically non-analytical behavior of the 
density field) is central, and that the simple context of 
non-expanding models should be sufficient to elucidate 
the essential physics. Note that we have not discussed 
here the amplitude of the correlation function, and in 
particular how it evolves in time, which is directly re- 
lated to the time evolution of the scale of non-linearity. 
This is where the fluctuations at large scales, which are 
different in the various IC considered, can play a role as 
envisaged in the cosmological context (through the linear 
amplification of power at large scales). We will address 
this question further, again considering non-expanding 
models, in future work. 
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